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ABSTRACT 



In this paper a new technique for modeling non-linear force-free fields directly from line of 
sight magnetogram observations is presented. The technique uses sequences of magnetograms 
directly as lower boundary conditions to drive the evolution of coronal magnetic fields between 
successive force-free equilibria over long periods of time. It is illustrated by applying it to MDI 
observations of a decaying active region, NO A A AR 8005. The active region is modeled during 
a 4 day period around its central meridian passage. Over this time, the dispersal of the active 
region is dominated by random motions due to small scale convective cells. Through studying 
the build up of magnetic energy in the model, it is found that such small scale motions may inject 
anywhere from 2.5 — 3 x 10^^ ergs s~^ of free magnetic energy into the coronal field. Most of this 
energy is stored within the center of the active region in the low corona, below 30 Mm. After 
4 days the build-up of free energy is 10% that of the corresponding potential field. This energy 
buildup, is sufficient to explain the radiative losses at coronal temperatures within the active 
region. Small scale convective motions therefore play an integral part in the energy balance of 
the corona. This new technique has wide ranging applications with the new high resolution, high 
cadence observations from the SDO:HMI and SDO:AIA instruments. 

Subject headings: magnetic fields - Sun:activity - Sunicorona 
Introduction 



The solar corona is a complex environment, 
where much of its complexity is due to magnetic 
fields. These magnetic fields are produced through 
a dynamo action near the base of the convection 
zone ( Charbonneau|2005 1 . Once formed they may 
become buoyantly unstable, rise through the con- 
vective zone and break through the photosphere 



(Archontis et al.|2004 


Magara|2004 


Galsgaard et 


al. 


2007 Murray & Hood 


2008 


Fan 


2009 


I where 



they are observed as either sunspots in white light 
or active regions in magnetograms. Such active re- 
gions structure the Sun's atmosphere and provide 
energy for eruptive phenomena such as solar fiares 



( |Bmi|[2008| and CMEs ( |Cremades eraL][2006 ). 

Once active regions form, small scale motions 
such as granular and supergranular flows result in 
the decay of the active region and the dispersal of 
magnetic flux across the solar surface in a random 
walk ( Leighton]|1964 ). This random walk leads to 
the convergence and cancellation of magnetic flux 
along Polarity Inversion Lines (PILs), in addition 
to the spreading of magnetic fields. Such evolution 
of magnetic elements in the photosphere acts as a 
driver for the build up of free magnetic energy in 
the solar corona. The subsequent evolution of the 
coronal field through quasi-static equilibria, in re- 
sponse to these motions, may be mathematically 
modeled through force-free magnetic fields, mag- 



1 



netic fields which satisfy j x B = where j 
( |Priest|[T982| . 



In recent years, non-hnear force-free field mod- 



eling has received much attention ( Schrijver et al 
2006[ [Metcalf et al.|[2008l ). A non-hnear force-free 



field, is a special class of force-free field, where 
the scalar function a{r) is a function of position, 
but must be constant along any field line. Non- 
linear force-free fields are of particular interest 



as they may contain free magnetic energy (Wolt 



jer 



1958). The construction of non-linear force- 
free fields (NLFFF) may be broadly split into 
two groups: static and time-dependent models. 
Examples of static modeling are the "extrapola- 
tion" of photospheric vector magnetic fields into 



the corona ( Schrijver et al.| 


2006 


Regnier & Priest 


2007 


De Rosa et al.| 


2009 


Wheatland & Regnier 


2009 


Jing et al.|2010 


, and direct fitting of NLFFF 



models to observed filaments and coronal struc 
tures (van Ballegooijen p004||Bobra et al.|2008[[&u| 



et al. 2009b |Savcheva fc van Ballegooijen 20091 



Static modeling involves constructing either indi- 
vidual NLFFFs or independent sequences of mag- 
netic configurations, where there is no direct cor- 
relation or evolution between the different fields. 

In contrast, time-dependent quasi-static mod- 
eling evolves the coronal magnetic field through 
continuous sequences of related non-linear force- 
free fields based on the evolution of a continuous 
time dependent lower boundary condition. This 
boundary condition may be specified through ei- 
ther idealized magnetic field configurations 



2003 



2009 



Ballegooijen et al. 2000 [Mackay fc Gaizauskas 



Mackay fc van Ballegooijen] [20051 |2006a 



or from observations (Mackay et al. 



2000 



Yeates et al!1|2007[ |2008[ ). Only through using dy 



namic models can the build up of free magnetic en- 
ergy and magnetic helicity in the corona be stud- 
ied, along with the evolution of isolated flux and 
current systems. Such a technique has been suc- 
cessfully applied in the past to model the evolution 



of the global corona ( 


Mackay & van Ballegooi- 


jen 


2006a|b 


Yeates & Mackay 


2009 


I, determine 



the origin and evolution of solar filaments (Yeates 



et al. 2008 Mackay & van Ballegooijen 20091, 



study the formation and lift off of magnetic flux 
ropes fMackay & van Ballegooijen 2006b Yeates 



origin of the Sun's Open Magnetic Flux ( 


Mackay 


& van Ballegooijen||2006b 


Yeates et al.||2010|. 



In this paper we apply time dependent quasi- 
static non-linear force-free modeling to model the 
effect that the dispersal of an active region has on 
the overlying coronal magnetic field. A key dif- 
ference between this paper and previous studies 
is the manner in which the time dependent lower 
boundary condition is applied. In previous stud- 
ies either idealized magnetic field distributions or 
simplified surface configurations derived from ob- 
servations were prescribed. In the present study, a 
new technique of prescribing the time dependent 
lower boundary condition is presented. The tech- 
nique allows the use of observed line of sight com- 
ponent magnetograms directly as lower boundary 
conditions, reproducing as observed, the disper- 
sal of the photospheric fiux of the active region. 
From this the effect of the surface motions on the 
coronal magnetic field and the subsequent energy 
input into the corona is determined. 

The active region chosen to illustrate this new 
technique is NOAA AR 8005. It is chosen as it 
was an isolated region and during the period of ob- 
servations there appears to be no significant fiux 
emergence and the fiux is well balanced through- 
out. The structure of the paper is as follows. In 
Section 2 the main properties of the decaying ac- 
tive region are discussed. In Section 3 the model- 
ing technique for both the photospheric and coro- 
nal fields are described. The results of the sim- 
ulation are given in Section 4, while in Section 5 
some simple calculations are compared to the sim- 
ulations to verify the results. Finally a discussion 
and conclusions are given in Section 6. 

2. Observations 

The decaying active region considered in this 
study is NOAA AR 8005. It emerged on the far 
side of the Sun and as it rotated onto the limb only 
dispersed magnetic polarities, without sunspots 
were present. In Figure [l] a full disk SOHO:MDI 
(Scherrer et al. 1995[ ) 96min line of sight com- 
ponent magnetogram from 18*'' Dec 1996 can be 
seen. The active region of interest lies at central 
meridian just below the equator. It has a simple 
magnetic morphology with a single positive and 
negative polarity. 

To study the evolution of the active region, 
SoHO:MDI 96min magnetogramf[^ are obtained 



^To produce the 96 min magnetograms, 5 individual 
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Fig. 1. — SoHO MDI 96 min line of sight compo- 
nent magnetogram for 18*'* Dec 1996 where white 
represents positive magnetic polarities and black 
negative magnetic polarities. The active region of 
interest, NOAA 8055, lies in the southern hemi- 
sphere at central meridian and is enclosed by the 
black box. 

from 19.12.05 UT on 16*'' December 1996 to 
20.47.05 UT on 20*'* December 1996, spanning 
a period of four days. The four day period was 
chosen to include, two days before and two days 
after central meridian passage, which occurs late 
on 18*'' December. During this period as LOS pro- 
jection effects are minimized, the measurements of 
the magnetic flux are known to be the most reli- 
able. In total 61, 96 min magnetograms cover the 
period of interest. The data are corrected for the 
area foreshortening that occurs away from central 
meridian using the IDL Solar Software routine 
drot-map. An area of 181 x 126 pixels is cut out 
of the rotated magnetograms centered on NOAA 
8055 where each pixel is 1.977". A time sequence 
of these can be seen in Figure [2] The area was 
chosen to be large enough to encompass the whole 
active region, but small enough that approximate 



SoHO/MDI magnetograms of higher time cadence and 
noise error of ±200 per pixel are averaged. Correspond- 
ingly the resulting 96 min magnetograms have a lower noise 
error per pixel of it9G. 



flux balance is achieved. The active regions that 
lie to the north and south of NOAA 8005 have no 
contribution to the flux in the area considered. 

On studying the magnetograms over the pe- 
riod of interest, the evolution of NOAA 8005 is 
dominated by the dispersal of magnetic flux from 
the main flux concentrations. There is no sig- 
nificant emergence of new magnetic bipoles. On 
studying Figure [2] or the movie of the evolution 
(moviel.mpg), the motion of the magnetic ele- 
ments is clearly dominated by a random walk due 
to supergranular flows and no strong shear or vor- 
tical flows are present. In Figure [3] the separation 
distance {S{t)) between the centers of flux in both 
the positive and negative polarities can be seen, 
given as a function of time (t). This separation 
is defined as S" =| S{t) \ where S{t) is the vec- 
tor pointing between the centers of flux in each 
polarity, 

(1) 

where Bz{i,j) is the line of sight component of 
the field at the i*'',j*'' pixel and R^j is the posi- 
tion vector of the i*'* J*'' pixel from the origin. The 
origin is defined to be the lower left corner of the 
magnetogram. Throughout the four day period, 
as the active region crosses central meridian, the 
separation of the polarities is seen to slightly in- 
crease from 68, 000km to 80, 000km. This increase 
is only 18% that of the initial separation. While 
in general, the region is diverging, there are sev- 
eral instances where patches of strong magnetic 
flux separate from the main polarities and cancel 
along the Polarity Inversion Line. A key point to 
note is that during the evolution no systematic 
shear motions or vortical motions of the magnetic 
flux are seen. The motions of magnetic elements 
are mainly due to the buffeting as a result of su- 
pergranular motions. Although supergranular mo- 
tions provide the only obvious flow pattern, due to 
the large latitudinal extent of the active region dif- 
ferential rotation may have a small non-negligible 
effect. The evidence and consequences of this will 
be discussed in Section 6. 

In Figure |4][a) the variation of the total un- 
signed magnetic flux within the selected active re- 
gion can be seen as a function of time. The flux 
values are given in units of lO^^Mx and time in 
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17-Dec-ig96 22:24-;05.423 ia-Dec-1S96 11:12:05.414- 




19-Dec-1996 01:39:05.405 19-Dcc-1S96 16:00:05.394 




-1996 04:51:05.3B5 



1996 19:11:05.375 



Fig. 2. — Time sequence of de-rotated SoHO MDI 
96 min line of sight component magnetograms of 
NOAA 8005. The time sequence covers the evolu- 
tion of the active region 2 days before and 2 days 
after central meridian passage (18*'' Dec 1996) 
where white represents positive magnetic polari- 
ties and black negative magnetic polarities. The 
images show the active region after it has been 
corrected for flux balance. 

days from 16*'* December 1996 19.12.05 UT. The 
solid line in Figure |4](a) shows the total unsigned 
flux as derived from the individual magnetograms. 
The flux values are consistent with a medium sized 
active region. Initially over the first day the flux 
values are seen to increase slightly, possibly due to 
LOS changes as the region rotates toward central 
meridian and the emergence of a small bipole on 
the south-western edge of the active region. How- 
ever, between days 1-4 the flux decays away as 
cancellation becomes signiflcant along the inter- 
nal PIL ([Green fc Kliem||2009[). By the end of 



day 4 the flux has dropped by 20% of its original 
value. This behavior flts well with the observa- 
tion that it is a decaying active region. The peak 
flux densities found within both the positive and 
negative polarities are around IkG. In Figure |4|^b) 
the flux difference (flux imbalance) calculated for 
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Fig. 3. — Separation between the center of flux of 
the positive and negative polarities as a function 
of time in days from 16*'' December 1996 19:12 
UT. 



each magnetogram can be seen as a function of 
time. Over the four day period the imbalance is 
systematically negative, possibly due to a domi- 
nantly negative surrounding small scale fields. At 
worst it is 10% that of the total flux of the active 
region. As this value is small, the active region 
may be regarded to be in flux balance to a high 
degree of accuracy. 

To use the magnetograms as a lower boundary 
condition in the numerical simulations complete 
flux balance is required. Therefore in each magne- 
togram the imbalance per pixel is determined and 
subtracted from each pixel. This correction varies 
from one magnetogram to the next. In general it is 
less than | 6 |G and so is comparable to the noise 
level within the magnetograms and signiflcantly 
less that the peak flux density. After correction, 
the variation of the total unsigned flux is given by 
the dotted line in Figure Qa) . It can be seen that 
this line follows the uncorrected flux very closely 
and applying the correction does not signiflcantly 
effect the net flux of the active region. 

3. The Model 

3.1. Coronal Model 

To simulate the evolution of the coronal mag- 
netic field above the active region a magnetofric- 
tional relaxation technique is applied. This tech- 
nique evolves the coronal magnetic field through 
a sequence of related non-linear force-free equilib- 
ria. The technique is similar to that described 
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Fig. 4.— (a) Total unsigned flux (lO^^Mx) as a 
function of time (days). The solid line denotes 
the value determined from the de-rotated data. 
The dashed line is once a correction has been ap- 
plied to ensure flux balance, (b) Imbalance of flux 
(10^°Mx) as a function of time (days). In both 
panels, time is from 16*'* December 1996 19:12 UT. 



in the papers of van Ballegooijen et al. (2000) 



and Mackay & van Ballegooijen ( 2006a I where the 



code described in these papers is adapted to a lo- 
cal cartesian frame of reference. Previously the 
technique has been successfully applied to con- 
sider the coronal field based on evolving photo- 
spheric boundary conditions in either idealized se- 



tups ( 


Mackay & Gaizauskas 


2003 


Mackay & van 


Ballegooijen|2005 


2009 


) or approximations to ob- 


served magnetograms ( 


Mackay et al.||2000 


Yeates 



et al.||2007[ 120081). A key element of these sim- 
ulations is that through considering a time se- 
quence of evolved equilibria, this allows the mem- 
ory of previous flux connectivities and currents to 
be maintained from one time to the next as energy 
and magnetic helicity is injected into the corona. 
Such a feature is significantly different compared 




Fig. 5. — Field line plot showing the initial po- 
tential field used in the simulation. The thin 
solid/dashed lines denote positive/negative con- 
tours of the normal field component on the base, 
while the thick lines are the field lines. 

to independent extrapolations which do not main- 
tain such memory. 

The magnetic field, B = V x A, is evolved by 
the induction equation. 



dA 

1h 



V X B 



(2) 



where v(r,i) is the plasma velocity and for the 
present study non-ideal terms are not included. 
To ensure that the coronal field evolves through 
a series of force-free states a magneto-frictional 
method is employed (Yang et al, 1986). We as- 
sume that the plasma velocity within the interior 
of the box (representing the solar corona) is given 

by, 



V = 



Ij xB 

B2 ' 



(3) 



where j = V x B and v is the coefficient of fric- 
tion. Equation ([s]) represents in an approximate 
manner, the fact that in the corona the Lorentz 
force is dominant and simulates the plasma expe- 
riencing a frictional force as it moves with respect 
to the reference frame. This ensures that as the 
coronal field is perturbed via boundary motions it 
evolves through a series of quasi-static non-linear 
force-free field distributions, satisfying j x B = 0. 
To carry out the computations a staggered grid is 
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used to obtain second order accuracy for all deriva- 
tives and closed boundary conditions are applied 
on the side and top boundaries to match those 
given in Section 3.2. 

A key new feature within this study is that 
the evolution of the coronal field is driven di- 
rectly through the insertion of a time sequence of 
observed line of sight component magnetograms. 
This means that the evolution of the photospheric 
field occurs in exactly the same way as is seen 
in the observations. No simplifications, idealiza- 
tion or smoothing is applied to the observed data 
and the data is inserted as observed (with only 
a small correction for fiux balance applied). The 
technique of specifying the sequences of the lower 
boundary conditions is discussed next. 

3.2. Lower Boundary Condition and Ini- 
tial Condition 

To model the evolution of the active region, 61 
magnetograms taken at discrete intervals of 96 
min are available, covering the four day period 
around central meridian passage. From these mag- 
netograms a continuous time sequence of lower 
boundary conditions is produced. This sequence 
of lower boundary conditions is designed to match 
each observed magnetogram, pixel by pixel, ev- 
ery 96 minutes. However, to model the evolution 
of the coronal field, A^t and Ay^ the horizontal 
components of the vector potential A on the base 
that correspond to this magnetogram must be de- 
termined. To determine these and produce a con- 
tinuous time sequence the following process is ap- 
plied, 

1. Each of the observed magnetograms, Bz{x, y, k) 
for /c = 1 — > 61 are taken, where k represents 
the discrete 96 min time index. 

2. Next the horizontal components of the vec- 
tor potential at the base, z = Q, are written 
in the form, 

A^b{x,y,k) = — , 
dy 

Ayb{x,y,k) = - — 
where $ is a scalar potential. 



3. For each discrete time index k, the equation 

^ _ dAyb _ dA^b 
^ dx dy 

then becomes, 

^ + ^ = _B (4) 

which is solved using a multigrid numeri- 
cal method. Details of this method can be 
found in the papers by Finn et aZ.(1994) and 
Longbottom (1998) and references therein. 
For a full description of the boundary condi- 
tions applied see |Mackay fc van Ballegooijen| 
( p009| . 

On solving for the scalar potential, $ this de- 
termines the horizontal components of the vector 
potential on the base {A^b, Ayb) for each discrete 
time interval, 96 min apart. To produce a contin- 
uous time sequence between each of the observed 
distributions a linear interpolation of A^b and Ayb 
between each time interval k and fc -I- 1 is car- 
ried out. To interpolate the fields 500 interpola- 
tion steps are used. By linearly interpolating the 
horizontal components of the vector potential on 
the base, this effectively evolves the magnetic field 
from one state to the other. Therefore techniques 
such as local correlation tracking are not required 
to determine the horizontal velocity. Numerically 
it also means that undesirable effects such as nu- 
merical overshoot or fiux pile up at cancellation 
sites do not occur and no additional numerical 
techniques to remove these problems have to be 
applied. 

The technique described above means that 
there are two time-scales involved in the evolution 
of the lower boundary condition. The first which 
is 96 min is the time scale between observations, 
the second which is 11.52s is the time-scale intro- 
duced to produce the advection of the magnetic 
polarities between observed states by interpola- 
tion along with the relaxation of the coronal field. 
The process described above reproduces the ob- 
served magnetograms at each 96 min discrete time 
interval and therefore produces a highly accurate 
description of the magnetogram observations. 

In movie2.mpg a comparison of the actual ob- 
servations (left) and normal magnetic field com- 
ponent produced by the above technique (right) 
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(a) 



(b) 



Fig. 6. — Field line plots showing (a) the non- linear force- free field after 4 days of evolution and (b) a 
potential field corresponding to the same normal field component as that of the non-linear force-free field 
after 4 days. 



can be seen. A good agreement between the two 
is found. The new technique means that neither 
feature tracking, nor Local Correlation Tracking 
is required to derive horizontal velocity profiles, in 
order to simulate the evolution of magnetic fields 
between observed distributions, therefore remov- 
ing a source of uncertainty. At this point it should 
be noted that the above technique only specifies 
Axb,Ayb on z = 0. It however does not specify 
Az which lies, 1/2 a grid point in the box and is 
determined by Equation [2j the coronal evolution 
equation. Non-potential effects near the base, as 
a result of the evolving lower boundary condition, 
may be contained within this term. These are sys- 
tematically produced as the horizontal field com- 
ponents on the lower boundary evolve in time. 

The choice and construction of the initial con- 
dition for the coronal field is now discussed. For 
the initial condition, a wide variety of choices may 
be made. These range from potential to non-linear 
force-free fields. One method to constrain this is 
to use coronal images of the active region to deter- 
mine the field geometry. The most realistic form 
of the initial condition is the non-linear force-free 
field. However, since no vector magnetic field data 
are available to us, such a field cannot be con- 
structed. On comparing various linear force-free 
field solutions with coronal images we find that no 
single value of the force-free parameter a fits the 



whole corona structure. Therefore as vector mag- 
netic field data are not available, we choose the 
initial condition to be a potential field, as the so- 
lution is both unique and stable. By choosing this 
initial condition, when analysing the evolution of 
the active region we restrict ourselves to studying 
trends in quantities, rather than absolute values. 
In future with SDOiHMI vector field observations 
this restriction of only considering trends will be 
removed. In Section 6 the consequences of chang- 
ing the initial condition from that of a potential 
field is discussed. 

To construct a potential field the equation 

V^A = 0, (5) 

is solved in a cube with sides ranging from 
< x,y,z < 6 on a 256'^ grid where 1 unit = 
60, 733km. The cube represents an isolated region 
of the solar surface, where flux may only enter or 
leave through the lower boundary (z = 0). As flux 
may only enter/leave through the lower boundary 
all field lines must start and end in the 2 = 
plane, and no field lines may leave through the 
side boundaries. This means that _B„, the normal 
component of B, vanishes on each of the faces ex- 
cept z — (which represents the photosphere). To 
satisfy this it is assumed that the tangential com- 
ponents of A are zero on each of the faces, except 
z = 0. In addition, the normal derivative of the 
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Fig. 7. — (a) Total magnetic energy (lO'^'^ Ergs) 
versus time from 16*'' December 1996 19:12 UT. 
The solid line represents the energy of a potential 
field, while the dashed line is for the non-linear 
force-free field, (b) Free magnetic energy (lO'^" 
Ergs) versus time. 

normal component of A is set to zero on each of 
these faces. If Equation ^ is solved subject to the 
above boundary conditions it is straight forward 
to show that the solution will have V.A = every- 
where within the domain (Finn et al., 1994). Thus 
the initial condition is a potential field associated 
with the imposed normal field on the boundary 
with the choice of the Coulomb gauge (V.A — 0). 
The initial potential field is constructed from the 
first observed magnetogram at 19.12.05 UT on 
16*'' December 1996. 

4. Simulation Results 
4.1. Magnetic Field Lines 

In Figure [5] an illustration of the field lines of 
the potential field used as the initial condition can 



be seen. From these field lines it is clear that 
the decaying active region has a simple bipolar 
form. The field lines connecting between the posi- 
tive flux (solid contours) and negative flux (dashed 
contours) form semi-circular loops. 

In Figure the illustration shows the field 
lines for the non-linear force-free field after 4 days 
of evolution. This figure can be compared to Fig- 
ure [6]jb) where the corresponding field lines for a 
potential field, deduced from the same normal field 
component on the base are plotted. In each case 
the starting point for the field lines is taken on the 
base within the positive flux region. On compar- 
ing the two images it is clear that there are many 
differences in the connectivity and structure of the 
fleld. The main differences occur low down along 
the PIL between the two main polarities. Here the 
fleld lines of the non-potential fleld have a much 
more sheared structure as a result of energy and 
helicity being injected along the fleld lines by the 
small scale convective motions. 

Another major difference is that for the fleld 
lines lying at the southern end, the connectivity of 
the field is very different from that of the poten- 
tial fleld. This is because the connectivity within 
the non-linear force-free simulation is initially de- 
flned at the start of the simulation and is preserved 
throughout the simulation (except where numeri- 
cal diffusion becomes large). The largest field lines 
within the simulation are only slightly different as 
the small scale convective motions have not been 
able to inject helicity along the full length of these 
field lines during the time period of the simulation. 

4.2. Magnetic Energy 

In Figure [Tj^a) the graph of total magnetic en- 
ergy stored within the coronal field can be seen as 
a function of time. The dotted line is for the non- 
linear force-free fleld simulation, while the solid 
line is for a potential fleld deduced from the same 
normal fleld component on the lower boundary as 
that of the non-linear force-free field (see Section 
3.2). Both values are initially equal to one-another 
as the initial condition is a potential field. For 
both cases the total energy is of the order of 10'^^ 
ergs which is typical for a small active region. 

Initially over the first day the energy of both 
the non-linear force-free field and potential field in- 
crease. This increase is due to the increasing flux 
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values, however after 1 day the magnetic energy 
in both peaks and then starts to decrease as the 
level of flux decreases due to cancellation. While 
both the non-linear force-free field and the poten- 
tial field have a similar behaviour, it can be seen 
that the non-linear force-free field has a system- 
atically higher energy compared to the potential 
field. This increased energy is due to the small 
scale convective motions which inject a Poynting 
flux into the corona and evolve the initial coronal 
fleld away from potential. 

In Figure a graph of the total free mag- 
netic energy within the coronal field is plotted as 
a function of time. The total free energy is defined 
as the difference between that of the non-potential 
and potential fields when integrated over the en- 
tire volume. This quantity is always positive and 
by the end of the simulation the free magnetic en- 
ergy is approximately 8 x 10'^° ergs which is 10% 
that of the potential field. An interesting feature 
is that throughout the simulation the free mag- 
netic energy as a result of the convective motions 
increases steadily. The rate of increase is around 
2.5 X 10^^ ergs s~^ indicating that the small scale 
motions deduced from the magnetograms may in- 
ject large amounts of energy into the corona. 

In Figure |8] the images illustrate the locations 
of free magnetic energy storage for days 1-4. The 
plots are in the x-z plane and give the free mag- 
netic energy summed along the line of sight. 



E{x,z) = A 



-dy 



(6) 



where B is the magnetic field of the non-linear 
force-free field and Bp the magnetic field of the po- 
tential field satisfying the same normal field com- 
ponents on the boundaries. The factor A repre- 
sents the area of the column being summed over 
{A = 2.02 X lO^*' cm^). The area factor is in- 
cluded so that the free energy along the line of 
sight is computed in units of ergs. In each plot the 
x-direction represents the full length of the com- 
putational box, but the z-direction is only half the 
height. White locations denote where the non- 
linear force-free field has a higher energy than the 
potential field when integrated along the line of 
sight (i.e locations where free magnetic energy is 
stored) . Black denotes where the non-linear force- 
free field has a lower energy when integrated along 
the line of sight (i.e locations where there is no free 



energy). We note that while the non- linear force- 
free field must and does have a volume integrated 
energy that is greater than that of the potential 
field (see Figure [t]), there is no restriction that 
in any sub-volume that this must always be true. 
Hence the negative values of the line integral quan- 
tity in Equation [6] are physically valid. 

As the small scale convective motions advect 
the field between days 1 and 4, the locations of 
free magnetic energy storage expand up into the 
coronal volume and mainly lie in the center of the 
box. From these images it is clear that the main lo- 
cations of free magnetic energy is low down in the 
corona (below 30 Mm) and between the two main 
flux concentrations. It is built up at this loca- 
tion, as here the magnetic fields are the strongest 
and the energy injection through Poynting Flux 
the greatest. In addition, these field lines are the 
shortest and occupy a smaller volume. 

While the images shown in Figure [8] are useful 
for illustrating the general spatial locations of free 
magnetic energy. Often such scaled images lead 
to misleading interpretation on the exact values. 
To quantify this, in Figure |9j graphs of the varia- 
tion of free magnetic energy can be seen for day 4, 
along the paths given by the dashed lines in Fig- 
ure |8] In Figure |9](a) the free magnetic energy can 
be seen as a function of x at a height of 7, 300 km 
{z = 0.12 units). Positive values denote locations 
where the non-linear force-free field has a higher 
energy, negative values where it is lower (when in- 
tegrated along the line of sight). From this graph 
it is clear that in the coronal volume the free en- 
ergy is mainly positive. Near the center of the do- 
main there is a small region where the non-linear 
force-free field has a lower energy than the poten- 
tial field. However the negative value of the line of 
sight integrated free energy is small compared to 
the locations where it is positive. From this it can 
be seen that the scaling and color shown in Fig- 
ure |8] exaggerates the level where the free magnetic 
energy is negative. In Figure [9]jb) the variation of 
free magnetic energy can be seen as a function of 
z, taken at a; = 3 units. It is clear that along this 
cut the free energy is always positive. In general 
the amount of free energy available decreases with 
height in the box. 

Alternative views of the free magnetic energy 



storage can be seen in Figure 10 In Figure [10[ a) 
it is shown in the y — z plane summed along the x- 
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Fig. 8. — Images showing the locations of free magnetic energy storage in the x-z plane summed along the 
y-direction for the non-linear force free field. White denotes locations where there is excess free magnetic 
energy compared to that of a potential field. The image is set to saturate at ±1 x 10^^ ergs. 



direction (Equation [7| and in Figure [lO[b) for the 
X — y plane summed along the z-direction (Equa- 
tion [s]). In each case results are shown for day 4 
and the same scaling is used as in Figure [Sj 



E{y,z)=A 
E{x,y)=A J 



(B2 - B 2 



8tt 
Stt 



(7) 



(8) 



Once again it is clear that the locations of free 
magnetic energy storage exist in the center of the 
computational box where the main flux concen- 
trations are at their strongest. In Figure [lOj^b) 
the main spatial location where the line of sight 
integrated free energy is negative is surrounded 
by positive free energy locations. To determine 



why such a region exits, Figure 11 compares the 
photospheric magnetic field distribution of the ac- 
tive region to the location where the line of sight 
free magnetic energy is negative. In this plot thin 
solid lines denote positive flux, dashed lines neg- 
ative flux and the thick solid line where the line 
of sight integrated free energy is —6 x 10^^ ergs 
or less. It can be clearly seen that the negative 
region is spatially correlated with the cancellation 
site between the positive and negative polarities. 
Therefore it is due to the removal of photospheric 
and coronal flux at this location as a result of flux 
cancellation (Figure |4]). 



4.3. Magnetic Helicity 

As the coronal magnetic field is subjected to 
the small scale random motions, these motions not 
only inject magnetic energy into the field but also 
magnetic helicity ( Demoulin fc Pariat|2009 ). This 
is seen by the increasing complexity of the field. 
To quantify this injection a calculation of relative 
magnetic helicity (Hr) is carried out. This quan- 
tity is defined as 



A.BdT 



Ap.Bpdr. 



(9) 



where A and B are the vector potential and mag- 
netic field of the non-linear force-free field. Simi- 
larly Ap and Bp are the vector potential and mag- 
netic field of a potential field which satisfies the 
same normal field components as that of the non- 
linear force-free field on all boundaries (namely 
Bz on z = and i3„ = on all other faces). The 
definition of relative magnetic helicity was first in- 
troduced by |Berger fc Field, ( ,1984, ) as an invariant 
and therefore meaningful measure of magnetic he- 
licity. It should be noted that in constructing the 
potential field for the relative helicity calculation 
the technique described in Section 3.2 is used. The 
horizontal components of the vector potential on 
the base, A^i, and Ayb, are identical for both the 
non-linear force-free and potential magnetic fields. 
This means that the expression for the relative he- 
licity does not include an addition surface integral 
term. 
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Fig. 9. — Graphs showing the free magnetic energy in units of 10^^ ergs as a function of (a) x at z = 0.12 
units and (b) z at a; = 3 units along the dashed hues shown in Figure [s] Positive values denote locations 
where the non-linear force-free field energy has a higher energy compared to that of a potential field. 



In Figure 12 the graph shows the variation of 
the relative magnetic helicity as a function of time 
(solid line). From the behaviour of the curve it 
can be seen that over the first day there is no 
definite trend of helicity injection by the surface 
motions. The relative helicity oscillates between 
positive and negative values. In contrast, between 
days 1-4 there is a definite trend of positive helic- 
ity injection, with the relative helicity showing an 
increasing positive value. Positive helicity injec- 
tion is consistent with observations that show that 
active regions in the southern hemisphere have a 



dominant positive helicity (Pevtsov et al. 19951 



The dash-dot line in Figurep^2|is a linear best fit to 
the relative helicity curve between days 1-4. The 
gradient of the line gives a helicity injection rate of 
1.218 X 10'^^ Mx^s~^. Such an increase of helicity 
within the coronal field can only be as the result 
of helicity injection through the lower boundary 
as a result of the relative motion of the magnetic 



fragments ( 


Demoulin & Berger|2003 


Demoulin & 


Pariat||2009 1 as during the period of observations 



there is no increase in magnetic flux. However, 
through considering the evolution of the magnetic 
elements in the magnetogram, there is no clear in- 
dication as to the origin of the helicity injection, as 
the dominant motion acting on the active region 
are small scale random motions. 



graphs of flux variation (Figure|4]) and separation 

between the centers of flux elements (Figure |3| , 
it can be seen that the positive helicity injection 
may possibly be related to large-scale properties 
of the evolution of the active region. To begin 
with on day 1 the flux peaks, after which it sys- 
tematically starts to decrease. This decrease is the 
result of first convergence and then cancellation of 
the polarities along the internal PIL of the bipole. 
Secondly, around the same time the separation 
of the sources increases, so there is an overall di- 
vergence of the active region polarities. Either of 
these two processes (convergence/divergence) may 
result in the net positive helicity injection. How- 
ever, within the present simulation we are unable 
to determine which. To compute this techniques 
such as Local Correlation Tracking would have to 



be applied ( Chae|2001 


Pariat et al.|2005 


Jeong & 


Chae 


20071. While there is no clear indication of 



In contrast, on comparing Figure 12 to the 



its origin, a third possibility is that such a positive 
injection of helicity is consistent with the effect of 
differential rotation as will be discussed in Section 
6. 

While the relative helicity has a positive value 
the amount of helicity is small. Active regions 
containing fiuxes of 1 x lO^^Mx may well contain 
helicity on the o rder of 10^^ _ iq^Mx'^ ( |Tian &: 
Alexander|2008 ) . On studying the coronal field in 
more detail it is found that the random motions 
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Fig. 10. — Images showing the locations of free 
magnetic energy storage on day 4 in (a) the y — z 
plane summed along the x-direction and (b) the 
X — y plane summed along the z-direction for the 
non-linear force free field. White denotes loca- 
tions where there is excess free magnetic energy 
compared to that of a potential field. The same 
scaling is used as in Figure [8| 

inject large amounts of both positive and negative 
helicity, but in near equal amounts over the four 
days of the simulation. This is what is expected 
from such motions, therefore resulting in a low net 
helicity. 



In Figure 13 a comparison of (a) the normal 



magnetic field component on the lower boundary 
and (b) the distribution of a on the lower bound- 
ary can be seen for the final time snapshot of 
the simulation. Similar plots were found at other 
times. In Figure [l^ a) white/black represents pos- 
itive/negative flux, while in Figure 13 'b) red/blue 
denotes positive/negative values of a. The mag- 
netic field values are set to saturate at ±100G 
while the a values saturate at ±2.73 x 10^^ m^^. 
On comparing the distributions it can be seen that 
the spatial distribution of a is well correlated to 
that of the magnetic field. The distribution of a 
take two forms. First, in locations of strong field 
concentrations the a distribution takes the form 



4.5 r 
4.0 r 
3.5 r , 
>. 3.0 r 
2.5 r 
2.0 r 
1 .5 r 




Fig. 11. — Comparison of the underlying mag- 
netic field distribution with locations were the 
non-linear force-free field has less energy than the 
potential field on day 4. In the plot positive and 
negative photospheric magnetic fiux is given by 
the thin solid and dashed contours. The thick solid 
contour denotes the area where the non-linear 
force-free field has an energy less that 6 x lO^^ergs 
of the corresponding potential field when summed 
in the z-direction. 

of large coherent patterns. Secondly, it takes the 
form of small isolated patches of intermingled sign 
which lie in small flux concentrations. It is clear 
from comparing the two images, that within the 
large concentrations of flux, a single magnetic el- 
ement of one polarity may contain both positive 
and negative values of a. Such a mixture of posi- 
tive and negative values of a, as a result of small 
scale convective motions, support the results that 
the coronal field of the AR cannot be modeled by a 
linear force- free field as found in Section 2.2. The 
positive and negative a distributions are similar 
to the divided distribution of current used in the 
idealised model of Regnier ( 2009 1 . An interesting 



feature of the comparison is that within the large 
flux concentrations the northern portions appear 
to contain more positive values of a, while the 
southern portions have more negative values. In 
addition to the a values that lie within the cen- 
ter of the simulation, values may be seen along 
the edge of the magnetogram area. Such edge val- 
ues are a pure boundary effect and should not be 
considered further. 

For the active region considered here no vector 
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Fig. 12. — Variation of relative magnetic helicity 
within the sinmlation as a fimction of time (sohd 
line). The dash-dot line denotes a straight line 
fit to the positive trend seen between days 1-4. 
The gradient of the straight line gives an average 
helicity injection rate over days 1-4 of 1.218 x 10'^^ 
Mx^s-^ 

magnetogram data are available in order to com- 
pare observed a distributions or trends of a in 
the observations with those obtained within the 
simulation. In future with the launch of SDO and 
the SDO:HMI instrument, regular vector magnetic 
field data will be available which can be com- 
pared to future simulations applying this tech- 
nique. Through this we will be able to deter- 
mine how much helicity is injected due to small 
scale motions, and in addition how much must 
be injected by other effects (e.g. torsional Alfven 
waves). This provides a new tool for the quan- 
tification of such effects which may be applied di- 
rectly to observations. 

5. Simple Calculations 

To determine if the rate of energy build-up 
within the present simulation is realistic, a sim- 
ple order of magnitude calculation is carried out 
to determine the rate of inflow of energy into the 
lower boundary via the Poynting Vector, S = 
(c/47r)E X B. The Poynting flux (P) through the 
lower boundary is. 




where E is the electric field and da = dan. In 
order of magnitude terms this equation may be 



simplified to, 

47r 

where Vh is the typical horizontal velocity, Bz a 
typical normal field component, Bh a typical hori- 
zontal field strength and A^n the area of energy 
input. The area of input is taken as the area 
of the observed magnetogram = 4.7 x 10^° 
cm^, Vh = O.^km.s^^ the peak flow rate of a su- 
pergranular cell. The values of the average field 
components are then taken to be B^ = 12G and 
Bfi = 2G, where these values are determined from 
the initial potential field. Through this an upper 
estimate of the Poynting flux becomes 4.5 x 10^^ 
ergs s~^. This value is around a factor of 2 greater 
than the value deduced within the simulations of 
2.5 X 10^^ ergs s~^. Although the order of mag- 
nitude calculation is larger by a factor of 2 this 
value is an over estimate as it assumes (i) the flow 
is always at the peak rate and (2) flows that are 
steady and systematic. In reality all elements in 
the magnetogram do not move at all times and in 
the simulations the flow is irregular. Taking this 
into account it can be seen that the simulation 
produces a realistic energy input into the coronal 
volume. 

In the second calculation, we consider how im- 
portant the energy input from the small convective 
motions is relative to an estimate of the radiative 
losses for the active region. The radiative losses in 
an optically thin plasma at coronal temperatures 
is given by, 

Er = nlR{T) ergscm-^s"^ (12) 

where rie is the electron number density, T the 
temperature and R{T) the radiative loss function 
which is given by Rosner et al (1978) as, 

R{T) w 10-17-73^-2/3^ g 3 < log T < 7.0 (13) 

The average temperature of the active region is 
determined to be log T = 6.41 through using the 
SXT filter ratio technique on the Al.l and AlMg 
filters. We take the volume of the coronal emitting 
plasma to be 10^^ cm'^. Note that this is not the 
volume of the computation box but rather an up- 
per estimate of the volume of emitting plasma in 
the active region. Taking also an upper estimate of 
the electron number density to be lO^cm"^, pro- 
vides an upper estimate of the coronal radiative 
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losses as ^ 1 X 10^^ ergs s~^. This value can be 
seen to be of the same order as the energy build 
up due to the small scale convective motions. 

While the two values are consistent with one- 
another, it should be noted that the free magnetic 
energy in the simulation is held within the coro- 
nal field and therefore not immediately accessi- 
ble to heat the corona and account for the radia- 
tive losses. For such a process to take place an 
additional energy release mechanism, would have 
to act. Possible mechanisms include nanoflares 



(|Parker|1988||Browmng et al.|2008j Wilmot-Smith 



et al. 20101 or turbulent reconnection(Heyvaerts 
fc Priest||1984 1. At this point it is unclear if such 
a mechanism may release the energy. Also the 
simple calculation does not take into account the 
effect of thermal conduction which is an impor- 
tant factor in the energy balance of hot coronal 
loops. However, taking these limitations into ac- 
count, the calculation does show that convective 
motions may play a key role in the energy build 
up of the corona. 

Using the rates of energy injection determined 
in the simulation and assuming that it remains 
constant over the life-time of the active region, it 
is possible to estimate how much free energy the 
active region could have on the 16th December, 
the start day of the present simulations, due to 
prior evolution as it rotated onto the disk and to- 
wards central meridian. Assuming that it emerged 
just before disk passage the active region would 
have evolved for 5 days. As a result the small 
scale convective motions would inject a total of 
1.1 X lO^^ergs of free magnetic energy into the coro- 
nal field. This value is 10% of the active region en- 
ergy on the 16th December. While the energy may 
have been 10% higher, we cannot use this to con- 
strain the initial condition without having vector 
magnetic field data to constrain the distribution 
of a on the lower boundary. 

6. Discussion and Conclusions 

In this paper a new technique for modeling non- 
linear force-free fields directly from line of sight 
component magnetogram observations has been 
presented. The key new feature of this technique 
is that sequences of line of sight magnetograms can 
be directly used as lower boundary conditions to 
drive the evolution of coronal magnetic fields be- 



tween successive force-free equilibria. As a result, 
the lower boundary condition in the model closely 
follows the observed magnetic field evolution as 
given by the observed line-of-sight magnetograms 
over the period of observations. 

The technique is illustrated by applying it to 
MDI observations of a decaying active region, 
NOAA AR 8005. The active region was close to 
central meridian and in flux balance to a high de- 
gree. It was modeled over a 4 day period, where 
61 individual SoHOrMDI 96 min magnetograms 
where available to use as lower boundary condi- 
tions. During this time the dispersal of the active 
region was mainly dominated by random motions 
due to small scale convective cells and flux cancel- 
lation along the internal polarity inversion line. 

Through applying the 61 magnetograms as 
lower boundary conditions and simulating a con- 
tinuous evolution of the surface and coronal field 
between them, the build up of free magnetic en- 
ergy into the corona was studied. It should be 
noted here that the initial condition for the coro- 
nal field in the simulation was taken to be a po- 
tential field. Due to the availability of only line-of- 
sight component magnetograms this was the only 
unique choice. However at the time corresponding 
to the start of the simulation the actual coronal 
field of the active region need not have been poten- 
tial. Therefore our aim was not to reproduce the 
evolution of the coronal loops as seen in coronal 
observations, but rather to study the trends of en- 
ergy input over the period of four days, where this 
is deduced from using observed magnetograms. 

When the initial coronal field is taken to be 
potential, the small-scale motions inject around 
2.5 X 10^^ ergs of free magnetic energy, this 
energy is mainly stored in the low corona below 
30 Mm. However, even if the initial state on the 
Sun was different, this would not significantly al- 
ter the trends of energy injection found within the 
simulation. To illustrate this we have repeated 
the simulation using an initial condition of a lin- 
ear force- free field, with a values of 1.55x lO^^m^^ 
and 1.86 x 10~®m~^. These values are 80% and 
88% of the first resonant value of a. Only posi- 
tive values of a are chosen as the active region lies 
in the southern hemisphere. With these alterna- 
tive initial conditions similar results and trends are 
found to those shown in Figure [7] The only differ- 
ence is that the magnetic field has a higher initial 
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Fig. 13. — Comparison of (a) magnetic flux distribution and (b) distribution of the force-free parameter a 
on the base (z — 0) at the end of day four. The magnetic field values are set to saturate at ±100G while the 
values of alpha saturate at ±2.73 x 10^^ m~^. 



energy. The subsequent rates of energy injection 
are 2.65 x 10^^ ergs s^^ and 2.8 x 10^^ ergs s~^ 
respectively. These values are consistent to those 
found when using the potential field. The slightly 
higher values occur due to the increased horizon- 
tal component of field in the linear-force free field 
leading to a higher Poynting flux. For all cases 
the energy injection through this new technique is 
consistent with estimates derived from poynting 
flux calculations. After just 4 days of evolution 
the free energy is 10% that of the potential field. 
Such energy buildup, is sufficient to explain the ra- 
diative losses at coronal temperatures within the 
active region (1 x 10^^ ergs s~^). Therefore small 
scale convective motions play an integral part in 
the energy balance of the corona. 

Through studying the evolution of the relative 
magnetic helicity in the coronal field, it was found 
that near equal amounts of positive and negative 
helicity are injected into the coronal field, but with 
a slight preference for positive helicity injection in 
the later stages. This is consistent with the ac- 
tive region being buffeted by small scale random 
motions and observations that show that active 
regions in the southern hemisphere have a pref- 



erence for positive helicity (Pevtsov et al. 19951 



One possible cause of this slight positive helicity 
injection is the effect of differential rotation. In 
the southern/northern hemisphere differential ro- 
tation will inject positive/negative helicity into the 
corona (iDeVorel 120001). iDemouhn et. all (|2002|) 



showed that helicity injection by differential ro- 
tation is composed of two terms, which for a de- 



cayed and dispersed active region are of opposite 
sign. The injection rate also increases with active 
region latitude and varies with active region tilt 
angle. During the short time period of our study, 
where the active region lies near the equator, we 
expect differential rotation to play only a minor 
role in ejecting positive helicity into the active re- 
gion. 

The new technique presented here has wide 
ranging applications for modeling the evolution 
of photospheric magnetic fields observed on the 
Sun and the subsequent effect these have on the 
coronal field over long periods of time. In future 
we shall use new high resolution, high cadence full 
disk vector magnetic field data from SDO:HMI 
in combination with coronal observations to con- 
strain the initial condition. With this we may 
then compare the evolution of the coronal mag- 
netic field with that seen in coronal observations. 
Under such constraints it will become a powerful 
tool to study the nature and properties of coronal 
magnetic fields. 
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